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The  Arbitrary  Ouasi-Ort.hogonal  Surface  Method  for 
Computing  Three-Dimensional  Flow  in  Turbine  Machinery* 


II.  Calculation  of  the  Three-Dimensional  Flow 
With  the  S, -surface  Twisted 

X 

Xin  Xiao-Kang  Jiang  Jin-liang  Zhu  Shi-can 
(Fudan  University) 


Abstract 


On  the  basis  of  the  auasi-orthogonal 
surface  method,  this  paper  obtains  a  three- 
dimensional  flow  field  solution  of  inviscous 
fluids  in  turbine  machinery.  This  method, 
under  the  condition  allowed  by  the  storage  of 
the  computer,  can  calculate  any  shapes  of  S, 
and  S2  surfaces  (including  twisting  - 

surfaces)  and  obtain  the  shapes  of  space  stream¬ 
lines  and  the  three-dimensional  distribution 
of  parameters  such  as  velocity.  The  calcula¬ 
tions  for  the  two  conditions  of  a  centrifugal 
compressor  with  a  small  impeller  (0110)  indi¬ 
cate  the  relative  value  of  the  twist  of  the  S,- 
surface  Is  about  3  -  10%.  The  S^-surface  x 

shows  a  double  twist  like  a  twist  drill  from 
the  inlet  to  the  outlet. 


Under  the  assumptions  that  the  fluid  is  incompressible, 
inviscous,  and  absolutely  non-rotational ,  Reference  [2]  used 
a  binary  solution  superposition  method  to  calculate  the  results 
with  twisted  and  S2  surfaces.  But  because  the  turbine  was 
not  real,  the  results  can  only  be  used  as  a  qualitative  refer¬ 
ence.  On  the  basis  of  quasi-orthogonal  surface  method,  this 
paper  eliminates  the  assumption  that  the  S1  surface  is  a  revolv¬ 
ing  body.  In  the  processes  of  integration  and  the  interpolation 
V - 
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by  equal  di  ffVrorvp  of  flux,  an  enui-narameter  transformation 
Is  introduced  to  transform  the  flow  region  on  the  equi-ortho- 
gonal  surfaces  to  a  rectangular  region  on  the  (f>  j?)  plane. 

It  is  then  possible  to  use  the  procedures  listed  in  Reference 
[1]  and  various  modification  procedures  to  obtain  the  iterative 
solution  of  the  entire  three-dimensional  flow  field.  In  the 
search  for  the  reverse  transformation  of  the  equi-parameter 
transformation,  this  paper  simultaneously  offers  three  numeri¬ 
cal  methods  (see  Section  III).  The  numerical  calculation  indi¬ 
cates  that  the  accuracy  is  very  good  in  the  reverse  transforma¬ 
tion  using  a  straight  line  quadrangle  instead  of  the  curved 
quadrangle.  It  is  an  approximation  method  which  drastically 
reduces  the  calculation  time.  Using  this  method  a  computer 
program  has  been  compiled  for  a  centrifugal  turbine  and  a  cal¬ 
culation  has  been  carried  out  on  the  modified  centrifugal  com¬ 
pressor  turbine  of  the  GJ11#  pressure  riser  of  Shanghai  Diesel 
Engine  Factory.  The  calculated  results  indicate  the  twist  of 
the  surface  is  relatively  small.  Therefore,  it  frequently 
uses  the  S1  surface  to  carry  out  the  calculation  for  the  revolv¬ 
ing  body  assumption  as  a  fairly  good  approximation.  The 
surface,  in  addition,  shows  a  double  twists  like  a  twist  drill 
from  the  inlet  to  the  outlet.  This  method,  in  principle,  is 
applicable  to  any  type  of  turbines.  However,  for  axial  flow 
machinery  we  must  add  inlet  and  outlet  corrections.  The  work 
in  this  area  is  currently  in  progress. 

I.  Basic  Assumption  and  Solution  to  the 
Velocity  Gradient  Equation 

This  paper  uses  tw.o  basic  assumptions:  1.  The  working 
medium  (or  fluid)  isan  inviscous,  equi-Cp,  and  equi-Cv  perfect 
gas.  The  flow  may  be  unisoentropic .  However,  for  ease  of 
description,  the  entropy  gradient  term  is  neglected  from  the 
equation.  The  friction  loss  of  the  real  gases  is  reflected  by 
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the  variation  process  same  as  the  one  used  in  Reference  [1]. 

2.  The  flow  passing  through  the  turbine  is  a  constant  relative 
to  the  turbine.  The  kinetic  equations  are  written  in  a  rela¬ 
tive  cylindrical  coordinate  system  . 


The  basic  equations  used  in  this  work  still  are  equations 
(1)  and  (2)  in  Reference  [1].  Here  we  slightly  modify  these 
equations  and  rewrite  them  as 


dw 

dq 


-  «  A II'  +  li  +  --  (--£  -  a,  *L\ 
I  W  \  dq  dq  / 


(1) 


where 


A 
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,  (  dr 

ds  \ 
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+  sin  a  sin  /?  cos  j? 


dO 

dq 


■  S-  ds  \  adW-  ,  ■  a  dr 

sin  a  H - cosa  Icospr - 3  —  2u>  sing  - 

dq  .  ,  .  dq  /  _  dm  dq 


„  (dWt  ,  ,  .  \  dO 

+  rcosp  I - *■  +  2m  sina  ) - 

\  dm  /  dq 


(2) 


On  the  quasi-orthogonal  surfaces,  two  independent  varia¬ 
bles  q  and  0  are  introduced  as  in  Reference  [1].  Therefore, 
m 

in  the  quqsi-orthogonal  surface  of  the  single  blade  channel  on 
the  (0,  q  )  plane  (it  is  the  region  between  the  pressure  sur¬ 
face  of  one  blade  to  the  suction  surface  of  another  blade)  the 
flow  field  is  shown  approximately  in  Figure  1.  The  intersection 
between  S2  surface  and  the  quasi-orthogonal  surface  is  indicated 
as  the  dotted  line.  The  intersection  of  the  surface  and  the 
quasi-orthogonal  surface  is  shown  as  the  point-solid  line. 

If  the  problem  we  are  solving  has  N  +  1  surfaces  and  H  +  1 
S2  surfaces,  then  the  equations  of  the  intersection  between  the 
quasi-orthogonal  surfaces  and  the  S2  surfaces  can  be  written  as 


0 -/«(*-)  (4-  1,2,  •••,«+  l) 


(3) 


The  interject. ion  equations  of  the  quasi-orthogonnl  rurfacor. 
and  the  S.^  surfaces  can  be  written  as 

'  ‘1-  “  0  -»  1,  2,  •••,  N  ■+■  1)  (i|) 

Obviously,  when  C  © )  =  const,  it  is  the  situation  that 
surface  is  a  revolving  surface. 

The  first  few  steps  of  the  problem  solving  procedure  used 
in  this  paper  are  similar  to  those  used  in  Reference  [1].  Steps 
1-5  are  basically  identical  to  those  in  Reference  [1].  The 
only  thing  is  that  the  coordinates  (r  ,  0  ,  z  )  of  an  arbitrary 
lattice  point  p  (it  is  the  intersecting  point  between  S.^  sur¬ 
faces,  S2  surfaces  and  the  quasi-orthogonal  surfaces)  on  the 
streamline  should  be  a  three-dimensional  numerical  set  and  it 
is  no  longer  a  two  dimensional  numerical  set  as  in  Reference 
[1].  Similarly,  a,  fi,  rt  ,  etc.  are  also  presently  a  three- 
dimensional  numerical  set. 


Figure  1.  Key:  1.  Pressure,  2.  Suction  surface, 

4.  Rotating  cover. 

With  respect  to  setp  6  in  Reference  [1],  we  carry  out  the 
following  modification.  In  order  to  make  point  p  or  the  stream¬ 
line  satisfy  equation  (1)  simultaneously  in  both  the  qm  and  the 
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0  directions,  the  following  modification  is  made  in  solving 
for  the  velocity  gradient  equation  as  proposed  in  Reference 
[1].  Assuming  that  the  value  W  at  point  A  in  Figure  2  is 
known,  then  we  first  proceed  with  a  numerical  liner  integration 
along  AB  to  find  out  the  W  value  of  all  the  points  along  AB. 

We  further  perform  integration  along  the  intersecting  line 
with  the  surface  to  obtain  the  W  value  at  all  the  lattice 
points  (Figure  2a).  We  shall  designate  them  as  the  VJ  [j,k] 
numerical  set.  On  the  other  hand,  let  us  start  from  point  A 
to  integrate  linearly  along  AC  in  order  to  obtain  the  W  value 
of  all  points  along  AC  and  then  integrate  along  the  inter¬ 
secting  line  with  the  surface.  The  obtained  W  values  at 
all  thelattice  points  are  designted  as  the  W^tjjk]  numerical 
set.  Finally,  let 
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<,)  <"> 

Figure  2 

us  take 

W[j.  it]  —  lililiiA]  +  //  “  1 »  2,  + 

2  •  H- t»2,  •••,//+ 1/ 

as  the  relative  velocity  values  of  all  the  lattice  points  on 
the  quasi-orthogonal  surfaces.  The  mathematical  meaning  of 
this  process  is  to  make  any  point  P  simultaneously  satisfy 
equation  (1)  in  two  directions  under  the  "average1'  idea  described 
above.  In  addition,  we  also  use  the  numerical  method  which 
starts  from  the  center  of  the  flow  field  toward  the  perimeter 
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to  r-'-'l-:  for  In  this  method  the  original  start  inf; 

toir.l  A  is  moved  to  the  center  of  the  flow  field  E  (Figure  3). 
This  would  further  raise  the  accuracy  of  the  numerical  solution. 


Figure  3 

II.  Flux  Integration  and  Equi-Parameter  Transformation 

The  major  difficulty  brought  about  by  the  elimination  of 
the  assumption  that  the  surface  is  a  revolving  surface  is 
how  to  carry  out  steps  8  and  9  used  In  Reference  [1],  From 
Reference  [1],  it  is  known  that  the  mass  flux  equation  in  the 
channel  between  two  neighboring  blade sthrough  any  quasi-ortho- 
gonal  surface  is 


.  c  “*  PgWmt o*(<A  —  a)rdOdqm  (5  ) 

where  the  integration  region  5  is  the  ABCD  region  on  the  quasi- 
orthogonal  surface  (Figure  2).  Since  the  surface  is  not  a 
solving  surface,  the  integral  of  equation  (5)  can  no  longer  be 
written  as  a  quadratic  integration  (i.e.  integration  first 
with  respect  to  0  and  then  to  q  }  as  in  Reference  [1].  It 
must  be  written  as  a  double  integral  as  in  equation  (5). 
However,  there  is  certain  difficulty  in  the  direct  integration 
as  in  equation  (5).  For  that,  we  have  to  introduce  an  equi- 
parameter  transformation  which  transforms  the  ABCD  region  on 
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the  surface  into  a  rectangular  region  on  a 

certain  transformation  plane  -  (if,?;)  plane  and  also  simul¬ 

taneously  transform  every  curved  quadrangles  regions  surrounded 
by  the  intersections  of  S.  and  S2  surfaces  into  rectangular 
region.  Thus,  the  regions  on  the  (;»  ;;)  plane  become  very 
regular.  The  intersection  line  of  the  surface  is  trans¬ 
formed  into  the  n  =  const  straight  line  and  the  intersection 
line  of  the  surface  is  transformed  into  the  £  =  const 
straight  line.  Therefore,  the  integrations  in  equation  (5) 
become  much  easier  on  the  Cl. »;)  plane. 

The  following  is  the  derivation  of  the  equi-parameter 
transformation  equations.  Let  us  assume  that  there  is  a  region 
ABCD  surrounded  by  a  curved  quadrangle  (Figure  4a)  and  the 
equations  for  the  four  sides  are  assumed  as: 

AB:0~U C^).  CI>:0 ■—/*;,(</«)  J  <6> 

We  are  seeking  for  a  transformation  which  transforms  the 
curved  quadrangle  ABCD  into  the  rectangle  A*B*C*D*  on  the 
(£>»;)  plane  (Figure  4b)  and  make 

A*C*:n  =  i),,  1 

W:S-S<,  (7) 

We  found  that  the  equi-parameter  transformation  below  can 
achieve  the  above  objectives: 


>c  ~ 

—  ■»« 
v  —  h  m 
—  V, 


/«+.G/~>  -  /*(</-)  ’ 

<*•<  »♦•<*> 


(8) 
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For  convenience 


e  can  simplify  equation  (0)  as 


Figure  4 

It  is  easy  to  write  the  Jacobi  determinant  which  transforms 
equation  (9)  as 
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Here  the  value  of  the  Jacobi  equation  /J/  at  the  lattice  point 
. <)  is 

1  7  ( /* .  >(  V-->.  O  —  0i .  *  I  f  it- .  t)_—  7 ■ .  *  1 

IJI’  1 1  -  /&/-,”«)  •  i-XJ.*)!  STi±,,'~"~  ( 13 ) 

where  Ai  -  £,„  -  ^  -  I///,  A>,  -  ijl+l  -  ,,  -=  l/.v. 

Using  this  quadratic'  numerical  calculation,  we  can  obtain 
the  calculated  flux  value.  In  order  not  to  treat  the  coordi¬ 
nates  "specially",  from  now  on  in  all  our  procedures  we  will 
use  the  "average"  method  as  we  did  in  solving  the  velocity 
gradient  ^quati^n .  F^r  example,  in  determining  the  flux,  we 
let 


G<  =  (  Jjl  f  7)  7-77  <*? 

Jo  Jo  111 

(14) 

(‘jf  CK(.S,  7)777  J'/ 

Jo  Jo  \J\ 

(15) 

then  the  calculate  flux 

G  —  (G,  +  C,)/2 

(16) 

Compare  G  with  the  given  fixed  flux  Gq,  if  \G  —  C0|/G,  <  ea  , 
then  switch  to  step  8.  Otherwise  use  the  method  similar  to  the 
one  used  in  Reference  [1]  to  correct  the  velocity  value  at 
point  A  or  E  and  repeat  steps  6  and  7. 

III.  Flux  Interpolation  by  Equal  Difference 
and  the  Reverse  Transformation 

For  step  8,  which  is  the  step  using  the  interpolation  by 
equal  difference  method  to  obtain  the  new  coordinate  positions 
of  the  streamline,  we  carry  it  out  on  the  Cf,  j?)  plane.  Because 
the  integral  is  the  quadratic  integral  form,  therefore  we  can 
proceed  using  the  same  method  as  in  Reference  [1].  First  we 
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obtain  a  set  ol’  new  streamline  coordinates  (£<;i»  17..1),. (  from 
equation  (4)  and  then  we  obtain  another  set  of  new  streamline 
coordinates  (.Li,  v.,)  from  equation  (15).  After  taking  their 
average  value,  we  get  the  new  positions  of  the  streamline 
coordinate  points  from  the  interpolation  by  equal  difference 
in  flux 


fo  —  (fC1  +  Sci)/2,  !Jc  —  +  lfcj)/2 


(17) 


The  next  problem  is  how  to  find  0^ 
physical  plane  (from  qmG  we  can  get  rQ 
streamline  coordinate  positions  £G  and 


and  a  „  on  the 
-mG 

and  zQ )  once  the  new 
ri  are  found . 


From  the  equi-parameter  transformation  (8),  we  know  that 
it  is  pretty  difficult  to  find  8q  and  qmQ  from  and  Pq.  It 
is  because  the  real  analytical  solutions  of  gXO )  and  /,(./„) 
are  actually  not  known.  We  have  used  various  methods  to  look 
for  the  reverse  transformation.  Here  we  are  only  going  to  dis¬ 
cuss  and  compare  three  methods  which  we  have  used  to  conduct 
detailed  calculations. 


1.  Rigorous  and  "accurate"  Reverse  Transformation. 

From  the  equi-parameter  transformation  (8),  we  know  that  when 
^  it;  ^  and  Vi  ^ 


£  *  „  _ d  —  )  . 

Vc.  ~  Vi  qm  -  gXO )  [ 

,,c  “  v, +,  -  y,  “  gl+X<»  -  g,(o)  1  ( 18 ) 


Our  purpose  is  to  locate  the  0Q  and  qmQ  which  satisfy 
equation  (18).  For  that,  we  are  rewriting  equation  (18)  as 


®  ^?»)  Sclf^+XqJ)  /*(</«,)]  1 

-  g, id)  +  yd  g,+XO)  -  g,W  I  I  (19) 
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If  an  initial  value  is  given  by  an  approximation 

method,  then  we  can  use  the  following  iteration  method  to 

obtain  the  iterative  values  of  9  and  q  : 

m 

“  lkU:">  +  UU .,(</!,•')  -  /*(-/!:'))  » 

9*-*+u  -  g,wml)  +  r  (20) 

Iteration  stops  when  I </!"*'’ —  </»>"’ I  <  .  Thus  the  reverse 

transformation  which  is  "accurate”  in  the  sense  of  approxima¬ 
tion  in  numerical  value  is  obtained. 

2.  The  Reverse  Transformation  Using  the  Linear  Quad- 
Rangle  Approximation. 

Because  it  is  time  consuming  to  use  the  above  equations,  the 
reverse  transformation  which  replaces  the  curved  quadrangle 
by  the  linear  quadrangle  is  proposed.  It  is  obvious  that, 
with  increasing  N  and  H,  this  approximation  is  very  good. 


Let  us  consider  an  arbitrary  curved  quadrangle  ABCD 
(Figure  4)  with  the  vertices 

u  *+l»  9*1+ ••*+>) 

The  reverse  transformation  equations  using  the  linear 
quadrangle  approximation  can  be  obtained  through  very  simple 
mathematical  operations  as 


e 

<r. 


bj  —  a  +  (c|  +  Ku,  0(f)?  —  d ) 
l-(rf + 

_  —  d  +  (fr)  +  K, —  a) 

I  —  (<-|  +  Kii.*X/q  +  K„t) 


-  flCf,  u) 
ImQ  ,  O 


(21) 


where 
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a  “  A’i,.  1 7"!  ■  i  ft.  4  •  ^  A||.  rfml.  *  ATli,  *4l7~l.  *41 

f  “  ^IM  ■+■  A'|,t  *41,  ^  A,.  *ft,  *  7„,.  ^ 

f  "t*  ^1*  Atl  /  ™  Aj4I.  4  ■*“  A|,  * 

~  ft.  «n  ft.  *  •  ^7«t  “  7»i+i.  <  qmi.  * 

#V„  ,  —  7 ■;  ' ■■  ♦  «_L  jaJ,  £  —  gg-'tl-  *t»  ~ 

^,.(*1  ft.*  ftn.*+ 1  $141.* 

»•  ^ .  «  '•  ♦  ft.  ♦  ^141.  »4I  ft.  *41 

“*  „  ^  —  ,1  ’  '•*■“  “  ,,  _ ,, 

V".I4I.  «  ‘/Ml.  |  7i»I  +  M4I  ¥•;.  *41 


Thus,  after  CG  and  nQ  are  known,  we  can  use  equation  (21)  to 

directlv  obtain  9-  and  q  _ .  The  results  of  numerical  calcula- 

“  u  mu 

tion  (shown  later)  indicate  that  this  method  almost  has  the 
same  precision  as  the  "accurate"  reverse  transformation.  How¬ 
ever,  the  calculation  time  can  be  drastically  reduced.  It  Is 
even  possible  to  reduce  the  computation  time  by  1/3. 

3.  Simplified  Uni-directional  Sampling  Approximation 

In  our  earlier  computation  process,  we  have  used  another 
simple  approximation  method  which  considers  the  reverse  trans¬ 
formation  approximately  as 

9»«(f),  7.  -  7«(*l)  (24) 


Therefore,  we  can  use  the  uni-direction  sample  strip  inter¬ 
polation  method  to  obtain  the  values  of  ft;*.*  and  ?»*;.*.«  from 
(£«,.*,  Vc.i.0  .  This  approximation  method  is  very  simple  but  its 
accuracy  Is  relatively  poor. 

After  the  reverse  transformation  is  obtained,  we  can  pro¬ 
ceed  with  the  same  steps  which  follow  as  in  Reference  [ID. 

The  iteration  continues  with  the  maximum  momentum  variation 
max  at  the  new  streamline  point  is  less  than  the  given  accuracy 

e  .  The  actual  details  of  program  compilation  are  omitted  here 
max  “ 


IV.  <~!nl  f'ulTl.f'ri  Examples  and  Analysis 


We  have  used  the  above  method  to  compile  the  719  computer 
program  of  our  school  and  carried  out  an  numerical  calculation 
on  the  modified  centrifugal  compressor  turbine  of  the  GJ11# 
pressure  riser  of  Shanghai  Diesel  Engine  Factory.  The  calcula¬ 
tion  regime  is:  inlet  temperature  =  293°K,  inlet  specific 
gravity  Y,  i.n.yskji/m'  ,  isobaric  specific  heat  Cr  —  1022. 33kg  • 
m/kg • °K  ,  variational  coefficient  nQ  =  1.506,  mass  flux  Cq  = 

0.3  kg/s,  angular  velocity  w  =  41 88.8  rad/s  ,  number  of  blades  = 
14,  number  of  quasi-orthogonal  surfaces  m  =  13^  number  of 
surfaces  N  +  1  =  9,  number  of  surfaces  H  +  1  =  9. 

Since  the  calculated  data  are  too  much  to  show,  we  are  only 
going  to  present  part  of  the  data  showing  the  major  relevant 
results.  Figure  5  shows  the  meridional  surface  lines  of  the 
turbine  and  the  meridional  line  of  the  S2  surfaces.  Figure  6 
gives  the  distribution  of  relative  Mach  number  on  the  average 
S2  surface  ( k  =  5)  and  also  shows  the  comparison  of  the  distri¬ 
bution  of  M  using  the  revolving  surface  assumption  with  our 
results  when  4  =  5,/“  1, 5,  y  (notice  that  the  m  coordinate  in 
the  figure,  one  uses  the  percentage  of  the  arc  length  of  the 
meridional  streamline  and  the  other  uses  the  percentage  of  the 
arc  length  of  spatial  streamline).  Figures  7-9  show  the 
distribution  of  M  and  the  diagrams  of  streamline  of  the  cross 
blade  spatial  streamlines  when  j  =  1,  5,  9.  Figure  10  indicates 
the  twisting  condition  of  the  intersection  between  the  quasi- 
orthogonal  surfaces  S1  and  S2  when  i  =  1,  5,  9,  13.  It  also 
simultaneously  gives  a  comparison  between  of  the  results  of 
calculation  obtained  using  the  revolving  surface  assumption  at 
the  inlet  surface  (  i  =  1)  and  those  using  the  method  described 
in  this  paper.  It  only  shows  the  deviation  of  the  surfaces. 
Figure  11  shows  the  comparison  of  results  obtained  using  the 
three  reverse  transformation  methods  on  the  calculation  of  the 
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twisting  of  the  Sn  surfaces  (only  1  =  1  Is  plotted).  Figure 
12  shows  the  condition  of  surface  twisting  when  the  number 
of  blades  Z.,  =  18. 

i  <i 

From  these  calculated  examples,  we  can  see  the  following: 

1.  From  Figures  6  -  9,  it  was  found  that  the  results 
obtained  using  the  method  in  this  paper  is  basically  similar 
to  these  obtained  using  the  method  described  in  Reference  [1]. 
From  Figures  10  and  12,  it  was  found  that  the  relative  twisting 
of  the  surfaces  is  relatively  small,  approximately  3  -  10". 
Therefore,  the  results  obtained  using  the  revolving  surface 
assumption  and  non-revolving  surface  method  are  very  close. 
Therefore,  in  general,  the  revolving  surface  assumption  is  still 
applicable . 

2.  From  Figure  10,  we  can  see  that  the  present  calculated 
shapes  of  the  S-^  surfaces  sometimes  appear  to  be  slightly 
convex  (the  surface  of  J  =  3  is  most  apparent).  The  height  at 
both  ends  varies  with  i  -  1,  5,  9,  13  differently.  At  the 
inlet  (i  =  1)  the  pressure  surface  q^  is  higher  than  the  suction 
surface.  This  phenomenon  is  consistent  with  the  results  obtained 
from  the  two-dimensional  simplification  of  the  axi-symmetrical 
flow  machinery.  At  the  front  middle  part  of  the  centrifugal 
turbine  (i  =  5),  the  S.^  surface  still  shows  the  phenomenon  of 
beinm  slightly  convex.  The  interseting  point  qm  on  the  pressure 
surface,  however,  is  lower  than  that  on  the  suction  surface. 
Therefore,  the  surface  shows  a  reverse  buckling  (twisting) 
phenomenon.  Further  back,  at  i  *  9,  13  (outlet),  another  reverse 
twisting  effect  appears  which  causes  that  the  qm  on  the  pressure 
surface  to  behigher  than  that  on  the  suction  surface  on  the  S.^ 
surface.  This  phenomenon  is  not  consistent  with  the  result 
obtained  in  the  two-dimensional  simplified  calculation  of  the 
axial  flow  machinery.  As  for  whether  this  belongs  to  the  special 
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Figure  10 
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Figure  12 

characteristics  of  the  centrifugal  machinery,  it  is  still  yet 
to  be  verified. 

3.  From  Figure  11,  it  was  found  that  the  linear  quadrangle 
approximation  seems  to  have  the  same  accuracy  as  the  "accurate" 
iterative  reverse  transformation  method.  The  relative  errors 
is  about  0.5%.  The  advantage  of  the  linear  quadrangle  approxi¬ 
mation  method  is  the  reduction  of  calculation  time  by  1/3.  For 
example,  when  the  convergence  accuracy  emax  =  1  mm,  both  methods 
require  16  iterations  while  the  calculation  time  can  be  reduced 
from  138  to  93  minutes.  Therefore,  the  use  of  the  linear  quad¬ 
rangle  approximation  is  reliable  and  time  saving.  The  use  of 
uni-directional  sampling  approximation  has  a  larger  error. 

The  maximum  relative  error  is  about  2%. 

Figure  12  shows  the  calculated  results  when  max  < 

-5 

10  ■'m.  At  this  time  the  maximum  momentum  variation  of  Mach 
number  is  3  x  10  .  This  can  be  considered  to  be  convergent. 

But  the  calculation  time  is  too  long. 
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V.  Conclusions 


This  paper  uses  an  equi-parameter  transformation  to 
further  develop  the  procedures  mentioned  in  Reference  [1]  in 
order  to  calculate  the  three  dimensional  flow  field  in  the 
blade  channel  for  an  inviscous  fluid.  The  condition  of 
twistin'  of  the  and  surfaces  can  also  be  quantitatively 
determined  simultaneously.  Based  on  the  results  from  one  tur¬ 
bine,  it  is  reliable.  But  due  to  the  lack  of  experimental  data 
and  analytic  solutions,  there  is  no  way  to  reach  certain  con¬ 
clusions  for  some  phenomena  and  results  (such  as  the  small 
de-tree  of  twistin'  of  the  surfaces,  the  double  twisting  pheno¬ 
menon,  etc.)  and  it  still  requires  further  computational  and 
experimental  verifications. 
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Abstract 

On  the  basis  of  the  author’s  previous  work1",  this  paper  deals  with  the  principles 
and  methods  for  the  combined  shaping  of  the  wing  root  regions  of  a  wing  body 
combination,  in  accordance  with  the  requirements  of  area  rule  and  shockless  pressure 
field. 

It  is  shown  that  the  favorable  non-midwing  arrangement  is  able  to  reduce 
the  amount  of  fuselage  alia  ng  required  for  the  reduction  of  the  zero  lift  wave  drag 
of  a  wing-body  combination.  The  combined  ahaping  >0  designed  may  thus  be  realized 
more  easily  in  structunl  arrangement. 


I.  INTRODUCTION 

In  fuselage  designs  based  on  the  "area  rule",  there  usually 
still  exists  the  problem  of  the  appearance  of  shock  waves  on  the 
wings  and  the  reduction  of  wave  drag  to  the  expected  value  cannot 
be  achieved.  If  one  requires  that  the  fuselage  shaping  be  such 
that  when  the  pressure  field  of  the  fuselage  is  added  to  that  of 
the  isolated  wing,  the  resultant  pressure  distribution  in  the  wing 
root  regions  will  be  congruent  with  the  pressure  distribution  on 
the  equivalent  yawed  wing,  then  shock  waves  will  definitely  not 
appear  on  the  wings.  In  this  manner,  the  "area  rule"  can  be  used 
to  full  advantage  [1-3].  The  fuselage  shaping  required  when  such 
a  pressure  field  method  is  adopted  is  generally  fairly  large.  In 
order  to  reduce  the  amount  of  shaping,  one  needs  to  resort  to  com¬ 
bined  shaping,  i.e.,  the  wing  shape  in  the  wing  root  regions  needs 
to  be  appropriately  tailored — the  wing  thickness  near  the  front  of 
1 - 
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the  v.’ir.j  rcet  r-..- -  Ions  has  to  be  Increased  and  the  position  of  max¬ 
imum  wing  thickness  has  to  be  shifted  gradually  from  outside  of 
the  wing  root  regions  towards  the  wing  root  regions  [4],  In  addi¬ 
tion,  we  have  found  that  if  a  favorable  non-midwing  arrangement 
is  adopted,  fuselage  shaping  can  be  substantially  reduced.  We 
have  given  in  reference  [5]  the  method  for  calculating  the  press¬ 
ure  field  of  the  fuselage  for  a  non-midwing  arrangement.  The 
equation  derived  there  will  be  directly  applied  in  the  present 
paper . 

2.  THE  AERODYNAMIC  DESIGN  OP  THE  WING  SECTION  AND  FILLET  IN  THE 
WING  ROOT  REGIONS 


For  a  given  Mach  number  M  in  supersonic  designs,  the  wing 
section  at  the  wing  root  for  a  given  pressure  distribution  can 
only  be  obtained  by  means  of  successive  approximations.  Such  a 
procedure  is  tedious.  The  key  to  achieving  faster  convergence  lies 
in  choosing  correctly  the  wing  section  in  the  wing  root  regions  in 
the  first  approximation. 


1.  Determination  of  wing  section  in  the  wing  root  regions 


in  the  first  approximation 


We  simplify  the  problem  by  considering  an  infinite  swept-back 
wing  of  constant  chord.  The  sweep  angle  A  is  equal  to  the  sweep- 
back  of  the  minimum  pressure  isobar  of  the  equivalent  yawed  wing. 
The  relation  between  the  linearized  pressure  coefficient  Cp  of  the 
infinite  swept-back  wing  under  the  conditions  of  subsonic  leading 
edge  and  the  slope  dz/dx  of  the  wing  section  has  been  given  in  [6J 
to  be 
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iijciv  cii'c  at  i'  1  u 1 t (_•  dx f f icuit its  associated  with  solving  for 
the  shape  of  the  win),’;  section  by  means  of  Equation  (1)  under  given 
pressure  dlr’ribuMon  on  the  equivalent  infinite  yawed  wing.  This 
Is  because  for  subsonic  critical  flow,  the  near  the  trailing 
edge  of  the  equivalent  infinite  yawed  wing  has  to  return  to  a 
positive  vain-.  .  .".vcord  i  nr  to  Equation  (1),  the  slope  of  the  wing 
section  should  be  positive.  This  causes  the  wing  section  to  stay 
open  along  the  trailing  edge.  Note  that  in  subsonic  flow  round 
infinite  swect-back  wing  of  constant  chord,  the  pressure  distribu¬ 
tion  on  the  wing  section  in  the  root  regions  is  given  by  [7,8]. 

II  d g  _  _ 

dx  1  ^  |"\/ 1  —  A/  -+■  tin  A  |  ds  J 

0*-*'  n  "l^/l  -  A/'cosM  -  sinAJ^J  (2) 

Comparing  Equations  (1)  and  (2),  one  sees  that  the  second 
term  in  Equation  (2)  is  nothing  but  the  pressure  coefficient  of  the 
wing  section  in  the  root  region  under  conditions  of  supersonic  flow 
round  an  object  (the  leading  edge  of  subsonic  flow).  Its  weight  in 
Equation  (2)  increases  with  the  M  number  of  the  oncoming  flow. 

Hence,  we  make  the  rressure  distribution  under  higher  subsonic  M 
numbers  (e.g.,  M  =  0.9)  equal  to  that  of  the  corresponding  equival¬ 
ent  infinite  yawed  wing  in  accordance  with  Equation  (2).  The  wing 
section  in  the  root  regions  thus  determined  not  only  results  in 
improved  pressure  distribution  on  the  wing  section  in  the  root 
region  under  transonic  flows,  but  also  causes  the  pressure  distri¬ 
bution  of  the  wing  section  in  the  root  region  to  converge  to  that 
of  Its  corresponding  infinite  yawed  wing  under  lower  supersonic 
flows.  Therefore,  it  is  reasonable  to  use  this  as  a  first  approx¬ 
imation  to  the  shape  of  the  wing  section  in  the  root  region  for  a 
given  M  number  in  supersonic  designs.  We  express  the  integral  and 
dz/dx  In  Equation  (2)  in  terms  of  sums  of  the  ordinate  z  of  the  wing 
section : 
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Thus,  for  a  riven  pressure  distribution  on  the  wing  root  wing 

section,  one  can  solve  from  Equation  (3)  the  ordinate  corres- 

Donding  to  any  specific  x  .  Here,  x  and  z  are  dimensionless 

U  U  U 

quantities.  The  characteristic  length  used  for  conversion  into 
dimensionless  quantities  is  the  wing  root  chord  c0-  In  this  set 
of  equations,  N  is  the  number  of  terms  in  the  sums.  Value  for 
*.»  4"  and  4‘J  corresponding  to  xy  in  the  cases  N  =  8,16,32  have 
been  tabulated  [7].  Suppose  the  variation  in  the  wing  root  wing 
section  will  converge  at  spanwise  position  •  Take  the  wing 
section  at  that  point  to  be  the  control  wing  section.  The  region 
between  the  wing  root  wing  section  and  the  control  wing  section 
is  called  the  wing  root  region.  The  spanwise  variation  of  the 
wing  section  in  the  wing  root  region  at  each  point  on  the  lines 
of  constant  percentage  chord  can  be  regarded  as  a  parabola  or  a 
cosine  curve.  As  the  entire  wing  region  within  the  front  Mach  line 
of  the  wing  section  in  the  root  regions  has  some  effect  on  the 
pressure  distribution  on  the  wing  section  in  the  root  regions,  one 
can  take  the  spanwise  coordinate  of  the  point  of  intersection  of 
the  said  front  Mach  line  and  the  leading  edge  of  the  wing  to  be 


Results  of  calculations  for  actual  examples  show  that  if  one 
makes  the  pressure  distribution  on  the  wing  root  coincide  exactly 
with  that  on  the  equivalent  yawed  wing,  negative  values  for  the 
ordinate  of  the  wing  root  wing  section  will  appear  near  the  trail¬ 
ing  edge.  This  means  that  the  upper  and  lower  surfaces  of  the 
wing  section  will  meet  at  some  point  in  front  of  the  trailing  edge 
and  not  close  to  the  trailing  edge.  This  is  not  reasonable.  Hence, 
one  can  only  match  the  pressure  distribution  on  the  wing  root  wing 
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Figure  1.  Chordwise  distribution  of  C  of  wing 
root  wing  section  P 

o  equivalent  yawed  wing 
x  original  wing  root  wing  section 

section  to  that  on  the  equivalent  yawed  wing  to  the  largest  poss¬ 
ible  extent.  Our  procedure  is  as  follows:  1)  Modify  the  press¬ 
ure  distribution  before  the  minimum  pressure  point  of  the  orig¬ 
inal  wing  root  wing  section  so  that  it  coincides  exactly  with 
that  on  the  equivalent  yawed  wing.  The  pressure  distribution 
after  the  minimum  pressure  point  should  be  between  the  pressure 
distributions  of  the  equivalent  yawed  wing  and  the  original  wing 
root  wing  section,  as  illustated  by  curves  P,  E,  F  in  Figure  1. 

2)  Make  the  pressure  distribution  before  and  after  the  minimum 
pressure  point  lie  between  those  of  the  equivalent  yawed  wing  and 
the  original  wing  root  wing  section,  as  illustrated  by  curves  H 
and  J  in  Figure  1.  Under  the  requirement  that  the  trailing  edge 
of  the  wing  section  should  be  closed,  the  closer  the  pressure  dis 
tribution  on  the  wing  root  wing  section  is  to  that  on  the  equi¬ 
valent  yawed  wing,  the  further  forward  will  the  point  of  maximum 
thickness  of  the  wing  root  wing  section  shift.  This  causes  the 
entire  wing  section  to  reduce  in  thickness.  To  overcome  this 
problem,  we  note  that  the  point  of  maximum  thickness  shifting  for 
ward  means  that  the  thickness  at  the  leading  edge  of  the  wing 


rccMcr.  in  *h-'  r-"  ~i onr.  needs  to  be  increased .  Hence,  v;e  can 

satisfy  this  requirement  without  having  to  r.ovc  the  point  of  max¬ 
imum  thickness  ef  section  excessively  forward  by  extending 

the  chord  of  the  wing  section  in  the  root  regions,  i.e.,  by  add¬ 
ing  a  fillet  to  the  wing  root  region. 

2 .  The  aerodynamic  design  of  the  wing  root  region  fillet 
in  the  first  aonroximation 


Consider  a  trapezoidal  swept-back  wing.  Suppse  the  position 
of  the  point  of  minimum  pressure  on  the  first  approximation  wing 
root  wing  section  has  already  been  shifted  from  point  C  on  the 
original  wing  root  wing  section  to  point  B.  Since  point  B  still 

falls  behind  the  point  of  minimum  pressure  A 
of  the  equivalent  yawed  wing  (Figure  2),  one 
needs  to  employ  the  method  of  lengthening  the 
chord  of  wing  root  to  make  the  point  of  minimum 
pressure  of  the  wing  root  wing  section  fall  on 
point  A.  Let  the  distance  between  point  B  and 
the  leading  edge  be  OB,  and  denote  its  ratio 
to  the  original  chord  OD  as  x^,  i.e.,  x-^  = 

Figure  2  OB/OD.  If  the  chord  of  the  wing  root  wing  sec- 

Key .  1  fillet  tion  is  extended  to  O'  such  that  the  point  of 

0 '  A  OB 

miminum  pressure  falls  on  point  A,  then  since  x^  =  and 

OA  =  OD*x  .  where  x  is  the  ratio  of  the  distance  OA  between  the 
o'  o 

point  of  minimum  pressure  of  the  equivalent  yawed  wing  and  the  lead¬ 
ing  edge  to  the  original  chord  OD,  one  obtains  for  the  extended 
chord  the  expression 


OD'  -  ’  (5 ) 

To  determine  the  spanwise  shape  of  the  leading  edge  of  the  wing 
root  region  fillet,  one  needs  to  determine  the  spanwise  variation 
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of  the  position  of  the  point  of  minimum  pressure  in  the  wing  root 
section.  Then,  in  accordance  with  the  requirement  that  the  line 
joining  the  points  of  minimum  pressure  on  all  the  wing  sections 
of  the  wing  root  regions  should  coincide  with  that  joining  the 
points  of  minimum  pressure  on  the  equivalent  yawed  wing,  the 
amount  of  chord  lengthening  required  for  each  wing  section  of  the 
wing  root  regions  can  be  determined  in  the  manner  described  above. 
Thus,  the  shape  of  the  leading  edge  of  the  fillet  in  the  first 
approximation  can  be  determined. 

Numerical  procedures  are  renerally  required  for  determining 
the  pressure  distribution  on  each  wing  section  in  the  wing  root 
regions.  For  the  design  of  the  shape  of  the  fillet  in  the  first 
approximation,  one  can  employ  the  method  suggested  in  C9J,  multiply 
with  the  first  term,  and  K2  with  the  second  term  on  the  right 
hand  side  of  Equation  (3)5  and  obtain 

c/x.) - — s  w*,  -  ( 6  > 
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In  the  above  equation,  A,.  is  the  sweep-back  of  the  line  joining 
points  of  constant  x^  in  the  wing  root  regions.  =  1  at  the 

wing  root  and  decreases  along  the  spanwise  direction.  For  its 
values,  see  Figure  1  of  [10].  For  the  values  of  K^,  see  Equations 
(10)  and  (11)  of  [10]. 

3.  ZERO-LIFT  SHAPING  FOR  THE  FUSELAGE  NEAR  THE  WING-FUSELAGE 
JUNCTION  FOR  A  GIVEN  M  NUMBER  IN  SUPERSONIC  DESIGNS 

After  the  calculations  in  the  first  approximation  have  been 
made,  the  pressure  distribution  and  fuselage  shaping  should  be 
computed  according  to  the  linearized  supersonic  theory  for  the 
three-dimensional  wing.  If  the  fuselage  shaping  Is  still  too 
excessive,  then  the  wing  section  in  the  wing  root  regions  obtained 
in  the  first  approximation  has  to  be  modified.  The  position  of 


maximum  Uiioxiitsi  of  ux-  wing  section  in  the  root  regions  should 
be  shifted  forward  to  beyond  15?  of  the  chord  and,  if  necessary, 
to  6.5-  of  the  chord.  The  ratio  of  the  ir.aximur.  thickness  of  the 
wing  section  to  the  chord  should  be  approximately  increased  so  as 
to  obtain  a  smooth  variation  along  the  chord  for  the  slope  of  the 
modified  new  wing  section.  The  pressure  distribution  and  fuselage 
shaping  should  be  again  computed.  The  process  is  repeated  until 
the  results  converge  to  what  is  desired. 

For  a  swept-back  wing-body  combination,  when  the  ratio  of 
the  wing  root  chord  to  the  radius  of  the  fuselage  is  greater  than 
6,  the  pressure  distribution  due  to  thickness  on  the  wing  of  the 
combination  is  basically  similar  to  those  on  the  isolated  wing 
obtained  by  regarding  the  side  of  the  fuselage  as  a  reflecting 
plane.  Of  course,  there  are  definite  differences  between  the  two 
in  the  wing  section  of  the  root  regions,  especially  in  the  rear 
part  of  the  root  chord.  But  this  will  not  cause  any  appreciable 
error  in  the  design  [1].  Hence,  one  can  regard  the  pressure  dis¬ 
tribution  due  to  thickness  on  the  wing  obtained  after  fuselage 
shaping  to  be  approximately  equal  to  the  sum  of  the  pressure  dis¬ 
tribution  on  the  wing  resulting  from  fuselage  shaping.  The  amount 
of  fuselage  shaping  is  determined  from  the  difference  between  the 
pressure  distribution  due  to  thickness  on  the  wing  and  that  on  the 
required  equivalent  yawed  wing.  Hence,  one  should  first  calculate 
the  pressure  distribution  due  to  thickness  on  the  isolated  wing. 
Then  calculate  the  pressure  distribution  resulting  from  the  fuse¬ 
lage  shaping  and  then  calculate  the  amount  of  fuselage  shaping 
required. 

1.  Calculation  of  pressure  distributions  due  to  thickness 
of  the  supersonic  isolated  wing 

Consider  a  thin  wing.  The  pressure  coefficient  due  to  thick¬ 
ness  at  the  singular  point  on  the  wing  plane  (z  =  0)  has  been  cal¬ 
culated  to  be  [11] 
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In  the  alcvv  equation  and  (ct,B)  denote  the  dimensionless 

quantities  corresponding  to  the  stationary  point  (x,y,l)  ana  the 
movin'-  point  (  f  ,  n  )  respectively: 

X  mm  f-~  JT*(  y  _  iL.  a  _  £  ~  g  mm  ~ 

After  this  coordinate  transformation,  the  front  Mach  line  equation 
£,E**±B(g  —  y)  for  point  (x,y,o)  on  the  z  =  0  plane  becomes 

.  -..,±0,,  [Xc(y)  +  x?(y)-x*(/,)j,  a,.miil^pt 

x,  de-nctes  the  dimensionless  quantity  corresoonding  to  the  coordi- 
1  *  _ 

nate  x^  =  x^/s;  c  denotes  the  dimensionless  quantity  corresponding 
to  the  wing  section  chord,  c  =  c/s;  s  denotes  the  half-span  of  the 
wing;  and  RqC^)  denote  the  radius  of  the  win?  leading  edge 

at  B-j-  and  62,  respectively;  Aq1  and  Aq2  denote  the  sweep-back  of 
the  wing  leading  edge  at  6-j^  and  62,  respectively,  where  B-^  is  the 
B  coordinate  of  the  point  of  intersection  of  the  right  front  Mach 
line  a  =  -  a2  of  the  fixed  point  (x,y,o)  and  the  right  leading 

edge,  and  32  is  the  3  coordinate  of  the  point  of  intersection  of 
the  left  front  Mach  line  a  =  +  a2  of  the  fixed  point  (x,y,o)  and 

the  left  leading  edge;  f(a,B)  is  the  wing  surface  slope  function  at 
point  (a,B),  V  M1-  1. 


For  a  subsonic  leading  edge  and  a  supersonic  trailing  edge, 
the  domain  of  integration  for  Equation  (7)  can  be  divided  into 
three  regions,  as  shown  in  Figure  3-  For  region  (1),  a  =  -  a^; 

for  regions  (2)  and  (5),  a  =  ai  +  a2 •  Thus  Equation  (7)  can 
be  rewritten  as 


Adad*~ii 

am,  -  j)  ] 

,  1  J  xJLtA  1 

*  ,  t  f  *,<M  1< 

b  +  tgA„ix  —  «*;(/?, )j 

(8) 


(a.  —  nV(fli  B )  t 

in  which  A  —  rr-r  ,,,7— - T - 77,*  and  .4  denotes  the 
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finite  portion  of  the  integral. 


As  the  integrand  A  in  Equation  (8)  contains  a  singularity  at 
a  =  0  (leading  edge)  and  another  one  at  a  =  +  a2  (on  the  front 

Mach  line),  these  singularities  have  to  be  removed  before  numerical 
integration  can  be  performed.  For  integrals  containing  square- 
root  singularities ,  such  as 
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we  can  make  the  following  variable  substitution: 

Let  a  —  o.25(a,  —  oj)(i  +  C)J(2  —  O.  .  The  integral  then  becomes 
(‘  in  which  the  integrand  is  positive  everywhere 

—  C,),/'  with  the  range  of  integration.  After  this  treat¬ 
ment,  no  more  singularities  exist  in  the  integrals  in  Equation  (8) 
and  the  Gaussian  method  can  be  applied  in  the  numerical  integra- 
t  ion . 
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Figure  3 


Figure  4 
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Let  the  equation  for  symmetrical  fuselage  shaping  be 


0K  v  '  ,  .  ,  „ 

•-  t  ■  t  )  CU5  2  tit) 

Ox  •  -  o 


(9) 


For  a  mid-win1''  arrangement,  the  pressure  distribution  on  the  wing 
and  the  fuselame  caused  by  symmetrical  fuselage  shaping  is  given 
by  [12,13].  (Let 
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(11) 


In  the  above  equations,  all  the  lengths  have  been  divided  by  the 
radius  Ro  of  the  cylindrical  section  of  the  unwaisted  fuselage; 
the  x-axis  lies  along  the  axis  of  the  fuselage,  the  y-axis  passes 
through  the  point  where  waisting  begins  and  (x,r,0)  form  a  cylin¬ 
drical  coordinate  system  (see  Figure  4);  R  is  the  radius  of  the 
fuselage  after  waisting;  t  is  a  specific  parameter,  usually  taken 
as  the  relative  thickness  of  the  wing  root  wing  section;  g2n  is 
the  amplitude  of  fuselage  shaping;  W2n  is  the  set  of  special  func¬ 
tions  associated  with  the  boundary  value  problem  of  solving  the 
wave  equation ;see  [14].  For  a  symmetrical  fuselage,  we  may  set 
n  *  0  in  Equations  CIO),  Cll).  Then  setting  r  =  1  will  give  the 
pressure  distribution  on  the  wing  root  wing  section  produced  by 
the  fuselage  shaping. 

Non-midwing  arrangement  is  an  Important  measure  for  reducing 
fuselage  shaping.  We  have  In  [5]  found  the  additional  pressure 


distribution  on  the  wing  root  wing  section  due  to  the  non-midwing 
arrangement  in  the  case  of  axially  symmetrical  fuselage  shaping 

tc  be 
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See  Equation  (39)  in  [5]  for  the  definition  of  the  functions  F  and 
N.  Hence,  the  pressure  distribution  produced  by  axially  symmetric 
fuselage  shaping  on  the  wing  root  wing  section  in  a  non-midwing 
arrangement  will  be  riven  by 
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Key;  1 — non-nidwing;  2 — midwing 
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3 •  The  design  of  fuselage  shaping  near  the  Junction  of 
the  wing  and  fuselage 


If  the  required  pressure  distribution  on  the  wing  root  wing 
section  is  the  wing  root  wing  section  pressure  distribu¬ 
tion  for  unwaisted  fuselage  is  and  the  pressure  distribu¬ 

tion  on  the  cylindrical  section  of  the  fuselage  due  to  the  effect 
of  the  forebody  is  C,,(*),  ,  then  we  require  that  the  pressure  dis¬ 
tribution  after  axially  symmetric  fuselage  shaping  be 
C,0O  —  Cfi(*)  —  C,,CO  —  For  a  non-midwing  arrangement,  this 

is  just  the  value  of  Cp  on  the  left  hand  side  of  Equation  (13) 
and  Equation  (13)  is  an  integral  equation  involving  gc(x)  which 
can  be  solved  by  numerical  integration  using  the  trapezoidal  for¬ 
mula.  If  we  take  the  interval  to  be  6,  then  when  x  =  m6,  by 
setting  j1  ..  —  <o(m  ~  la-)t  one  obtains  from  Equation 
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In  the  above  equation,  r  (0)  can  be  determined  from  Equation  (13) 
to  be 
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In  the  case  of  a  midwing  arrangement,  let  F  =  N  =  0  in  Equations 
(14)  and  (15).  After  obtaining  the  value  for  gQ(m6),  we  integrate 
Equation  (9)  for  the  condition  that  n  =  0  by  means  of  the  trap¬ 
ezoidal  formula  by  an  approximation  and  obtain  the  radius  of  the 
waisted  fuselage. 

K(mff)  —  1  +  I  ga({)^{ 
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If  we  require  that  the  fuselage  be  smooth  where  the  wing-fuse¬ 
lage  junction  commences,  i.e.,  that  f0(0)-\0,  then  we  see 

from  Equation  (15)  that  we  must  have  Cp(0)  =  0  at 

x  =  0.  If,  in  addition,  we  require  the  axial  cross-section  area 
distribution  satisfy  the  "area  rule",  then  the  fuselage  shaping 
will  not  be  axially  symmetric.  Assume  the  fuselage  shaping  equa¬ 


tion  to  be 
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and  the  area  distribution  demanded  by  "area  rule"  to  be  S(x). 
Then  Csee  Figure  5). 


From  Equation  (18)  one  obtains 
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Suppose  the  pressure  distribution  produced  by  fuselage  shaping 

as  required  by  the  wing  root  wing  section  is  t(C  (x):  then  from 

po 

Equations  (11)  and  (12),  we  obtain  (considering  that  AC^  is  approx¬ 
imately  also  applicable  in  the  axially  symmetrical  case) 
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Substituting  Equation  (19)  into  the  above  equation,  one  obtains 
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in  which  we  let 
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(22 ) 


The  integrals  in  Equations  (21)  and  (22)  can  be  computed  by 
approximate  integration  using  the  trapezoidal  formula.  After  sub¬ 
stituting  g  (.x)  from  Equation  (19)  and  determined  from  Equa¬ 

tion  (21)  into  Equation  (17),  one  can  obtain  the  amount  of  fuse¬ 
lage  shaping  AR(  x,0)  at  different  positions  x,8.  For  the  case 
of  a  r.idwing  arrangement,  let  N  =  F  =  0  in  each  of  the  above 
equations . 

Finally,  the  fuselage  shaping  at  the  wing- fuselage  junction 
must  be  such  that  it  smoothly  connects  to  the  fuselage  shaped 
according  to  the  "area  rule"  behind  the  trailing  edge  of  the  wing 
root  wing  section. 

The  results  of  computation  for  actual  examples  show  that  for 
a  given  reduction  in  zero-lift  drag,  the  amount  of  shaping  is 
smaller  for  a  lower  isolated  wing  arrangement  and  is  smaller  tnan 
that  required  for  a  midwing  arrangement.  This  is  easy  to  compre¬ 
hend  as,  in  the  non-midwing  arrangement,  the  shielding  effect  of 
the  fuselage  on  the  upper  surface  of  the  wing  is  enhanced,  thus  aug¬ 
menting  the  effect  of  the  fuselage  shaping.  Hence,  non-mldwing 
arrangement  is  an  important  measure  for  reducing  the  amount  of 
fuselage  shaping. 
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Abstract 

The  free  diffractions  o:  a  weak  shook  wave  issuing  trom  an  open  end  of  a  shock 
tube  and  the  formation  and  development  of  a  vortex  ring  from  the  end  of  the  tube 
are  observed  by  a  schlieren  method.  Also  observed  are  the  oblique  reflections  from 
the  horizontal  floor  and  the  vertical  plate. 

As  the  wave  described  above  impinges  on  the  horrizontal  surfuce  near  grazing 
incidence,  the  Mach  stem  forms.  Under  certain  conditions,  when  this  Mach  stem 
travells  along  a  slant  surface  with  an  incidence  angle  which  is  larger  than  the  critical 
value  with  respect  to  the  horrizontal  axis,  the  phenomena  of  "dissolution”  of  the 
Mach  stem  occurs. 

The  observation  of  such  a  new  phenomena  and  its  analysis  are  presented  in  this 
paper. 

In  1971,  we  used  a  Schlieren  method  to  study  th^  flow  dia¬ 
grams  of  a  weak  shock  wave  emitted  from  a  shock  tube  and  the 
induced  ultra-pressure  field.  The  shock  wave  inside  the  tube  is 
planar  and  the  gas  flow  velocity  in  the  wake  of  the  wave  is  in  the 
subsonic  range.  After  issuing  from  the  opening  of  the  shock  tube, 
the  wave  Itself  is  diffracted  and  undergoes  regular  reflections 
from  the  horizontal  floor  and  the  vertical  plate.  At  the  same 
time,  the  gas  flow  in  the  wake  of  the  wave  forms  a  vortex  ring 
outside  the  tube  which  travels  forward  at  a  definite  velocity.  In 
addition,  a  sonic  emitted  flow  field  Is  formed  near  the  open  end  of 
the  tube.  When  this  shock  wave  travels  along  a  slanted  surface, 
one  further  observes  the  so-called  "dissolution"  phenomenon  in 
which  the  incident  wave  and  the  reflected  wave  travel  separately 
along  the  slanted  surface.  Such  "dissolution"  is  a  new  phenomenon 
which  we  will  describe  and  explain  on  a  physical  basis  in  this 
paper.  Observations  have  been  made  on  the  flow  diagrams  of  waves 
issuing  from  both  porous  and  nonporous  tube  walls. 
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WAVL  DIAGRAMS 


After  the  weak  cheek  wave  issues  frer.  the  open  end  of  the 
shock  tube.  It  is  first  diffracted.  The  corner  of  the  open  end 
of  the  tube  can  be  rer-arded  as  a  source  of  interfering  waves. 

Hence,  the  shock  v;ave  will  be  attenuated  by  these  interfering 
waves,  and  its  wavefront  changes  from  a  planar  to  an  ellipsoidal 
shape  when  it  has  traveled  a  distance  equal  to  the  distance  of  the 
tube  (plate  I,  photo  1). 

The  main  shock  wave  undergoes  regular  oblique  reflection 
from  the  horizontal  floor.  Under  our  experimental  conditions,  one 
can  calculate  the  variation  of  the  intensity  of  the  shock  wave  with 
distance  from  the  tube  opening  by  measuring  the  angles  of  incidence 
and  reflection  (plate  I,  photo  2).  One  should  note  that,  although 
the  reflection  from  the  floor  is  an  oblique  reflection,  its  des-- 
tructive  effect  is  not  negligible.  Roughly  speaking,  for  medium 
strong  shock  waves,  the  ultra-pressure  in  the  wake  of  the  reflected 
wave  is  on  the  same  order  of  magnitude  as  that  of  the  main 
reflected  wave.  For  weak  shock  waves,  it  could  even  exceed  the 
latter . 

For  reflection  of  the  shock  wave  from  a  vertical  plate  placed 
at  a  distance  equal  to  twice  the  diameter  of  the  tube,  the  rela¬ 
tion  between  the  angles  of  incidence  and  reflection  is,  on  the 
whole,  that  for  regular  reflection  (plate  I,  photo  3).  After  the 
shock  wave  is  reflected,  it  suffers  from  the  disturbance  from  the 
vortex  ring  and  the  emitted  gas  flow  and  is  no  longer  a  spherical 
wave.  Rather,  it  becomes  bow-shaped  in  the  two-dimensional  case, 
and  saucer-shaped  in  the  three-dimensional  case.  After  it  passes 
through  the  vortex  ring,  the  main  reflected  wave  increases  the 
pressure  in  the  original  flow  region.  In  the  case  of  a  porous 
wall,  the  reflected  wave  is  still  saucer-shaped. 
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In  the  process  of  backward  movement  along  the  floor  of  the 
Incident  wave  reflected  off  the  vertical  plate  and  the  reflected 
wave,  Mach  reflection  Is  sure  to  occur.  The  height  of  the  Mach 
stem  of  this  Mach  reflection  (i.e.,  the  height  of  the  fork)  in¬ 
creases  as  the  wave  moves  on,  and  the  intensity  of  the  Mach  wave 
gradually  decreases.  Therefore,  under  the  conditions  of  given 
incident  wave  intensity  and  given  angle  of  incidence,  the  height 
of  the  Mach  stem  can  be  qualitatively  used  as  a  basis  for  estimat¬ 
ing  the  attenuation  of  the  Mach  wave.  Under  our  experimental  con¬ 
ditions,  the  position  of  this  fork  rises  with  a  slope  of  approx¬ 
imately  15°.  The  wave  diagram  for  reflection  off  the  vertical 
plate  is  shown  in  plate  I,  photos  4  and  5. 

For  the  case  where  the  vertical  plate  is  placed  at  a  distance 
equal  to  twice  or  four  times  the  tube  diameter,  the  reflection 
near  the  tube  is  a  Mach  reflection  whether  the  wall  is  porous  or 
not.  Experimental  results  show  that,  under  these  conditions,  the 
destructive  effects  of  the  reflected  wave  on  the  obstacles  and 
structures  in  the  vicinity  of  the  wall  is  obviously  larger  than 
that  of  the  diffracted  main  shock  wave. 

2.  THE  EMITTED  GAS  FLOW  AT  THE  OPENING  OF  THE  TUBE 

Two  cases  need  to  be  considered  for  the  emitted  gas  flow,  viz. 
the  case  where  the  flow  after  the  shock  wave  in  the  tube  is  sub¬ 
sonic  and  the  case  where  It  is  supersonic. 

For  the  case  where  the  gas  flow  at  the  cross-section  of  the 
tube  opening  is  Initially  subsonic,  generally  speaking,  before 
the  vortex  ring  leaves  the  end  of  the  tube,  the  divergent  waves 
arising  at  the  corner  of  the  tube  opening  will  be  transmitted  into 
the  tube.  This  results  in  increased  velocity  of  the  gas  flow  at 
the  cross-section  of  the  tube  opening,  which  eventually  becomes 
sonic.  The  wavefront  of  the  above  mentioned  divergent  waves  have 
a  stronger  thinning  effect . 


Th*'  o  m  i  °  i  *n.r  f!  fol  lowing  the  min  "heck  wave  gradually 
develops  from  the  unaffected  trumpet-shaped  region  into  the  peach- 
shaped  rvgi'U  as  a  result  of  the  thinning  effect  of  the  diverr.ent 
waves.  This  peach-shaped  region  is  where  the  two  divergent  waves 
arising  fin::,  the  corner  of  the  tube  opening  interfere.  It  should 
be  pointed  out  that  despite  the  fact  that  the  main  shock  wave  is 
continually  impinged  upon  by  the  divergent  waves,  as  the  gas  flow 
emitted  from  inside  the  tube  keeps  on  replenishing  energy,  the 
ultra-pressure  of  the  shock  wave  in  the  vicinity  of  the  central 
region  of  the  gas  flow  is  actually  larger  than  that  in  the  regions 
above  and  below.  This  is  characteristic  of  the  diffraction  of  the 
shock  wave  after  issuing  from  the  tube  opening. 

After  the  gas  flow  at  the  cross-section  of  the  tube  opening 

has  been  accelerated  to  sonic  velocity  by  the  divergent  waves,  as 

the  pressure  at  the  tube  opening  is  higher  than  that  outside  the 

tube,  the  gas  flow  will  continue  to  be  accelerated  beyond  the  tube 

opening,  forming  a  local  supersonic  flow  region  at  the  tube  opening 

(plate  I,  photo  6,  plate  II,  photo  7).  When  the  ultra-pressure  of 

2 

emission  is  1  kg/cm  ,  the  Mach  number  measured  for  the  local  super¬ 
sonic  flow  for  a  tube  with  a  nonporous  wall  is  1.18;  while  that  for 
a  tube  with  a  porous  wall  is  1.05.  Therefore,  the  Mach  number 
decreases  approximately  by  11?  when  the  wall  is  porous. 

For  the  case  where  the  gas  flow  at  the  cross-section  of  the 
tube  opening  is  supersonic,  supersonic  emitted  flow  is  formed  at  the 
corner  of  the  opening  and  the  divergent  waves  from  the  corner  will 
not  be  able  to  enter  the  tube.  As  the  pressure  at  the  tube  opening 
is  larger  than  the  atmospheric  pressure  of  the  surrounding,  the 
free  boundary  of  the  gas  flow  will  certainly  expand  outward.  The 
flow  diagram  is  similar  to  that  for  the  stable,  constant  supersonic 
emitted  flow. 
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VORTEX  RIRO 


Because  cf  the  el't'ect  of  adhesion,  when  the  gas  flows  around 
a  pointed  corner,  a  small  kernel  of  vortex  ring  with  a  diameter 
just  a  little  larger  than  the  inner  diameter  of  the  tube,  is 
first  ferry' d  close  to  the  end  of  the  tube.  As  time  roes  on,  the 
gas  behind  the  main  shock  wave  keeps  feeding  the  vortex  ring. 

The  vortex  ring  grows  fatter  and  fatter  until  it  finally  leaves 
the  tube  opening  and  travels  forward  at  a  definite  velocity. 
Because  the  vortex  ring  undergoes  turbular  flow  diffusion  and 
momentum  exchange  with  the  surrounding  gas  flow,  it  further 
increases  in  diameter  and  fatness.  However,  as  a  result,  the  vor¬ 
tex  ring  becomes  looser  and  looser  until  it  finally  falls  apart. 


After  the  vortex  ring  leaves  the  tube  end,  it  is  first  accel¬ 
erated.  The  velocity  reaches  a  maximum  when  it  has  traveled  a 
distance  approximately  twice  the  tube  diameter.  What  follows  is 
decelerated  motion,  shattering  and  dissipation. 


For  shock  waves  driven  by  high  Dressure  compressed  air  in  a 

2 

nonporous  tube,  when  the  ultra-pressure  of  emission  is  1  kg/cm  , 

the  average  velocity  of  the  vortex  ring  is  127  m/s;  when  the  ultra- 

2 

pressure  is  q  kg/cm  ,  the  average  velocity  of  the  vortex  ring  is 

156  m/s.  For  shock  waves  driven  by  explosive  gas,  when  the  ultra- 

2 

pressure  of  emission  is  15  kg/cm  ,  the  average  velocity  of  the  vor¬ 
tex  ring  is  129  m/s. 


For  shock  waves  driven  by  high  pressure  compressed  air  in  a 
porous  tube,  because  the  emitted  shock  wave  is  attenuated,  the 
velocity  of  the  vortex  ring  is  lower  by  about  10£.  In  the  case  of 
the  nonporous  tube,  the  time  elapsed  before  the  vortex  ring  leaves 
the  tube  opening  is  about  250  ys.  In  the  case  of  the  porous  tube, 
the  vortex  ring  leaves  the  tube  opening  earlier  by  about  8£. 
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4  . 

OF  THE  POKES 


WALL  ALD  THE  OISSIl ATIHG  EFFECT 


After  the  main  shock  wave  passes  through  the  pores  in  the 
wall,  it  is  emitted  out  of  the  pores  and  diffracted.  The  waves 
leaking  out  of  the  norer,  diffuse  rapidly  and  their  intensity  is 
rapidly  attenuated.  At  the  same  time,  these  waves  superpose  on 
one  another  and  form  a  Mach  cone  envelope  with  finite  amplitude 
(plate  II,  photo  8).  The  oblique  reflection  from  the  floor  and 
the  vertical  plate  of  the  main  diffracted  wave  has  similar  wave 
diagrams  for  porous  and  nonporous  walls  (plate  II,  photo  9).  The 
pore  content  of  the  porous  walls  discussed  in  this  paper  is  l.^h% 
or  2.7853.  The  pores  are  uniformly  distributed. 


What  comes  after  the  main  shock  wave  in  the  tube  is  a  high 
pressure  gas  flow.  Under  the  action  of  this  high  pressure  gas  flow 
the  gas  particles  will  certainly  be  emitted  through  these  small 
pores  into  the  stationary  gas  outside  of  the  tube.  Hence,  a  small 
emitted  flow  is  formed  at  each  pore.  When  the  shock  wave  is  barely 
emitted  from  the  tube  opening,  the  emitted  flow  field  near  the 
pores  has  already  formed.  Approximately  at  the  time  the  vortex 
ring  is  at  a  distance  equal  to  1.5  times  the  tube  diameter  from 
the  tube  opening,  owing  to  the  adhesion  effect  of  the  gas,  the 
turbulap^Jlow  diffusion  and  momentum  exchange  among  the  gas  parti¬ 
cles  in  the  center  region  of  the  emitted  flow  and  those  in  the  sur¬ 
rounding  stationary  gas  result  in  the  merging  of  the  small  emitted 
flow  fields  at  their  outer  edges. 


When  the  reflected  shock  wave  passes  through  the  emitted  flow 
region  at  the  tube  opening  and  the  emitted  flow  region  near  the 
pores,  the  shock  wave  will  develop  two  wavefronts  (.plate  II,  photo 
10). 


The  waves  escaping  from  the  pores  can  be  nothing  more  than 
waves  with  finite  amplitudes.  It  is  thus  completely  reasonable  to 
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deduce  the  intensity  ol'  emission  of  the  main  shock  wave  at  the 
tube  openin'-'  from  the  ingle  of  the  Mach  cone  of  the  envelope  of 
the  „aves  escaping  through  the  pores.  The  shock  wave  intensity 
obtained  this  way  agrees  fairly  well  with  that  measured  with 
ether  mo  the  Is. 

The  analysis  given  in  Section  1,  2  and  3  shows  that -the  pores 
in  the  tube  wall  can  indeed  reduce  the  intensity  of  the  emitted 
shock  wave.  Experimental  results  show  that  the  intensity  of  the 
emitted  shock  wave  is  less  for  porous  walls  than  for  nonporous 
walls  by  about  16$  (at  the  tube  opening)  and  9%  (at  a  distance 
1.3  times  the  tube  diameter  from  the  opening).  Qualitatively 
speaking,  we  consider  the  drop  to  be  about  10$. 

5.  OBSERVATION  OF  A  WEAK  SHOCK  WAVE  TRAVELING  ALONG  A  SLANTED 
SURFACE  USING  OPTICAL  METHODS  AND  ITS  ANALYSIS 

When  the  shock  wave  emitted  from  the  end  of  the  tube  and  the 
wave  reflected  off  the  floor  travel  together  along  a  slanted  sur¬ 
face,  the  flow  diagrams  and  the  formation  of  the  wave  system  is 
very  complicated.  First,  the  Incident  shock  wave  and  the  wave 
reflected  off  the  floor  are  not  planar  but  are  curved.  Next, 
before  they  hit  the  slanted  surface,  they  have  already  formed  Mach 
reflection  but  with  a  rather  short  Mach  stem.  Furthermore,  after 
they  hit  the  slanted  surface,  the  interfering  wave  arising  from 
the  corner  first  impinges  on  the  Mach  stem  and  the  original  re¬ 
flected  wave.  In  general,  after  the  Mach  stem  hits  the  slanted 
surface  at  a  definite  angle  of  incidence,  reflection  of  the  Mach 
stem  is  expected  to  occur.  However,  we  found  in  our  experiments 
that  in  the  actual  flows  the  Mach  stem  "dissolved”  and  reflection 
of  the  original  Incident  wave  and  that  of  the  wave  reflected  off 
the  floor  occurred  separately.  This  is  evident  in  many  of  the 
pictures  taken  (plate  III,  photos  13-20). 


Now  let  us  explain  this  phenomenon.  First,  we  analyze  the 
wave  diagrams  for  the  case  where  the  Mach  stem  has  been  formed 
before  the  weak  shock  wave  hits  the  slanted  surface.  This  is 
the  case  for  the  30°  and  4b 0  slants. 

Rirht  after  the  shock  wave  hits  the  corner  of  the  slanted 
surface,  the  disturbance  of  the  corner  results  in  the  formation 
of  a  series  of  compression  waves.  As  the  gas  flow  behind  the 
shock  wave  is  always  subsonic  relative  to  the  shock  wave,  the  com¬ 
pression  waves  will  keep  impinging  on  the  shock  wave.  The  wave- 
front  of  these  compression  waves  is  the  fan-shaped  arc  in  the 
photos.  As  the  shock  wave  is  very  weak,  the  reflected  wave  of  the 
Mach  stem  merges  with  the  wavefront  of  the  compression  waves. 

Note  that  the  formation  and  growth  of  the  Mach  stem  in  the  Mach 
reflection  is  always  caused  by  the  incident  wave  being  pushed  from 
the  back  by  the  oblique  reflected  wave .  In  our  case,  the  slope  of 
rise  of  this  Mach  stem  is  about  15°.  Hence,  the  process  of  travel¬ 
ing  along  the  30°  and  45°  slants  is  simply  the  process  of  reflection 
of  the  Mach  stem  from  the  slanted  surface  and  that  of  the  short¬ 
ening  of  the  Mach  stem.  Those  parts  of  the  entire  wave  (the  orig¬ 
inal  incident  and  reflected  waves)  that  have  not  been  disturbed 
or  reflected  by  the  slanted  surface  or  the  corner  will  continue 
traveling  in  the  original  direction  at  the  original  speed.  As  the 
waves  progress  along  the  slanted  surface,  the  Mach  stem  is  grad¬ 
ually  shortened  and  eventually  the  fork  meets  the  slanted  surface. 
What  follows  is  the  individual  reflection  and  shortening  of  Mach 
stem  of  the  incident  and  reflected  waves.  The  reflected  wave  of 
the  original  wave  is  tangent  to  the  wavefront  of  the  compression 
waves  Issuing  from  the  corner;  the  reflected  wave  of  the  original 
reflected  wave,  on  the  other  hand,  merges  with  the  wavefront  of  the 
compression  waves.  This  Is  how  the  Mach  stem  Is  "dissolved"  and  how 
the  original  incident  and  reflected  waves  separate  to  be  independ¬ 
ently  reflected  and  to  travel  separately  along  the  slanted  surface. 


[Plate  111,  photos  13-15  (.30°  slant);  Plate  III,  photos  16-18 
(^5°  slant)] 


For  the  wave  diagrams  in  which  the  Mach  stem  has  not  been 
formed  before  the  incident  shock  wave  hits  the  slanted  surface, 
the  process  of  traveling  along  the  slanted  surface  is  the  same  as 
the  foregoing,  i.e.,  one  in  which  the  incident  wave  and  the  re¬ 
flected  wave  are  separated  (plate  IV,  photos  19-20.  (60°  slant)). 
However,  for  a  15°  slant,  as  the  slope  of  the  slanted  surface  is 
about  the  same  as  that  of  the  rise  and  growth  of  the  Mach  stem, 
the  incident  wave  and  the  reflected  wave  do  not  separate  in  the 
process  of  traveling  along  the  slanted  surface.  Here  one  only 
observes  the  impingement  of  the  compression  waves  on  the  original 
shock  wave  (plate  II,  photos  11,12). 

Thus  we  see  that  the  wave  diagrams  for  the  weak  shock  wave 
traveling  along  a  slanted  surface  depends  on  the  Intensity  of  the 
wave,  the  slope  of  rise  of  the  Mach  stem  and  the  fork  and  also  on 
the  slope  of  the  slant  Itself.  To  obtain  the  detailed  quantitative 
relationship  among  these  factors,  one  needs  to  make  further  exper¬ 
imental  observations  and  perform  the  corresponding  numerical  cal¬ 
culations  . 
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